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We perform a statistical analysis of the binary black hole problem in the post-Newtonian approx- 
imation by systematically sampling and evolving the parameter space of initial configurations for 
quasi-circular inspirals. Through a principal component analysis of spin and orbital angular mo- 
mentum variables we systematically look for uncorrelated quantities and find three of them which 
are highly conserved in a statistical sense, both as functions of time and with respect to variations 
in initial spin orientations. We also look for and find the variables that account for the largest 
variations in the problem. We present binary black hole simulations of the full Einstein equations 
analyzing to what extent these results might carry over to the full theory in the inspiral and merger 
regimes. Among other applications these results should be useful both in semi-analytical and nu- 
merical building of templates of gravitational waves for gravitational wave detectors. 
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I. INTRODUCTION 



The inspiral and merger of binary black hole systems 
is of great importance in many astrophysical settings, 
from the large-scale structure of the universe via galactic 
merger trees to earth-based gravitational wave detectors 
(LIGO Q], Virgo 0, and GEO600 [3]) wherein stellar- 
mass black hole inspirals and intermediate-mass black 
hole mergers are the best source candidates in their fre- 
quency band. Because they encounter a small signal-to- 
noise ratio, these detectors rely on banks of templates of 
gravitational waves for their detection, which vary based 
on the properties of their sources. Within General Rela- 
tivity (GR) the parameter space of initial configurations 
for binary black holes in quasi-circular orbit is eight- 
dimensional. From the no-hair theorem each black hole is 
uniquely characterized by its mass rrii , spin orientation S^ 
and dimcnsionless spin magnitude ^ (z = 1, 2) defined by 
spin vector Sj = ^imfSj. In the absence of matter there 
is no preferred scale in GR so that the parameter space 
can actually be trivially reduced to at least seven dimen- 
sions. Still, the cost of generating templates for a multi- 
dimensional parameter space through full numerical sim- 
ulations is prohibitively expensive. Semi-analytical mod- 
els are a very promising approach but they also require 
calibration through numerical simulations, which are still 
computationally expensive even if to a lesser extent. Any 
possible simplification or guidance in building or comput- 
ing gravitational wave templates is therefore crucial. 

The post-Newtonian (PN) approximation to General 
Relativity is a good description of the inspiral dynam- 
ics until the black holes get close to each other [H |5J. 
Because the PN approximation is dramatically simpler 
than full GR — and correspondingly much more compu- 
tationally inexpensive — it is feasible to thoroughly study 
the parameter space in this approximation numerically 



and obtain information that is relevant in the fully rel- 
ativistic case. Among recent developments this includes 
statistical predictions for recoil velocities after a black 
hole collision [6 , how large recoils can be suppressed due 
to spin alignment with the orbital angular momentum 
[7], whether initially uniform spin orientation distribu- 
tions remain uniform during evolution, how initial and 
final spin orientations correlate over time, and the use 
of graphics processing units to accelerate these simula- 
tions [5]. 



In this paper we present an analytical and numerical 
statistical analysis of the binary black hole parameter 
space in an initially quasi-circular orbit in the PN ap- 
proximation. We numerically sample and evolve the en- 
tire 7-dimcnsional space of initial configurations in the 
PN approximation. We focus on finding the most and 
least relevant variables in the spin dynamics, namely, 
those that account for most of the dynamics and those 
that are largely conserved, even in the presence of radia- 
tive corrections in the binary's evolution. This is pre- 
cisely the area of Principal Component Analysis (PC A), 
a rather standard technique in multivariate statistical 
analysis that does not yet seem to have been applied 
to the binary black hole problem. To obtain approx- 
imate expressions in closed form we also implement a 
PCA in what we call the instantaneous approximation 
and validate to what extent the numerical results can be 
described via this analysis. Finally, we present some pre- 
liminary studies, through numerical simulations of collid- 
ing binary black holes using the full Einstein equations, 
aimed at studying to what extent the results of this pa- 
per carry over to the fully relativistic case in the inspiral 
to merger regimes. 



II. POST-NEWTONIAN APPROXIMATION 



averages over the unit sphere 



We work with the PN equations from Ref. (9J[T0], which 
describe a quasi-circular inspiral of two spinning black 
holes up to 3.5PN order in the angular frequency uj and 
spin effects up to 2PN order with the covariant spin sup- 
plementary condition. The evolution is given by a system 
of coupled ordinary differential equations (ODEs) for the 
orbital frequency w, the individual spin vectors Sj of each 
black hole and the unit orbital angular momentum vector 
L. For completeness, these equations are given explicitly 
in Appendix A. The spin vectors Sj are given in terms of 

their magnitude and orientation by Sj = rn^Xj^j where 
Xj € [0, 1] is the dimensionless Kerr spin parameter of 
the j th black hole. 

We numerically evolve the system of coupled ODEs for 
w, L and Si from some initial frequency uji to a final one 
uif and systematically sample the range of masses and 
spin parameters rrij , Xj ■ We typically choose u>i corre- 
sponding to an initial separation of r ~ 40Af and a final 
frequency of up to ujf = 0.05 (with that exact final fre- 
quency unless otherwise stated), which is a conservative 
estimate of where the PN equations still hold [H [5] . See 
[H] for details on our numerical implementation of these 
equations. 

In the PN approximation the free parameters of the 
problem have different hierarchies [HUH] because the spin 
orientations (Si, S2) evolve as a function of time while 
their masses and spin magnitudes (mj,Xj) are constant. 
For this reason it is appropriate to consider the spin 
orientations as stochastic variables for each rrij ,Xj [H] • 
Since there is no preferred scale in vacuum GR, in all 
of our numerical simulations we fix the total mass to 
M = mi + m 2 = 1. Also, we typically use, unless other- 
wise stated, 40,000 random spin orientations to capture 
patterns in evolutions for each tuple (toi , xi > m 2, X2)- We 
have found that this is sufficient to capture the main fea- 
tures of our results. 



III. PRINCIPAL COMPONENT ANALYSIS 

In PCA one seeks to determine the variables that are 
statistically relevant and to dimensionally reduce from 
the problem those that are not. The covariance between 
two stochastic variables X, Y is given by 



1 

-in 



(■■■)dtt, 



(3.2) 



Cov(X, Y) 



{(X - (X))(Y - (Y))) 
(XY)-(X)(Y) 



(3.1) 



and provides a measure of the degree to which their fluc- 
tuations are correlated. A smaller (larger) covariance im- 
plies lower (higher) correlation. In particular, the covari- 
ance of a variable with itself is its variance (i.e., the stan- 
dard deviation squared) and measures deviations from 
the mean value. The brackets above represent expecta- 
tion values. In the context of this paper they represent 



i.e., expectation values with respect to the two black hole 
spin orientations. 

When there are multiple stochastic variables Xi (i = 
l,...,n) one can construct their associated covariance 
matrix C with components Cij = Cov(Xi,Xj). This 
matrix is symmetric, non-negative definite and can be 
diagonalized with an orthogonal transformation. If we 
denote the components of the i th normalized eigenvector 
Vj by V? then the principal components (PCs) are the 
associated eigenmodes, 



£i = V>X X + ... + V t n X n . 



(3.3) 



The PCs are uncorrelated, a consequence of the orthogo- 
nality of the eigenvectors, and their associated eigenval- 
ues A are their variances, 



Cov{£i,£j) = XiSi. 



(3.4) 



The fact that, by construction, principal components are 
uncorrelated with each other is important since they pro- 
vide independent pieces of statistical information. 

The smaller an eigenvalue Xi the more likely that the 
corresponding linear combination £i will not deviate from 
its average value for a randomly chosen pair of spin ori- 
entations. Therefore, if there exist small eigenvalues then 
the associated PCs are largely conserved in a statistical 
sense. Conversely, the larger an eigenvalue is then the 
more relevant the associated PC is in describing the dy- 
namics and variations in the problem. 

There are two related but different senses in which a 
principal component with small variance will be shown to 
be semi-conserved in this paper. The first is in the usual 
sense, i.e., being constant as a function of time for an 
arbitrary but fixed initial spin configuration. The second 
is in the statistical sense that deviations of a principal 
component from the mean value are small for arbitrary 
but fixed initial and final times over a set of runs with 
the same masses and spin magnitudes but random spin 
orientations. In other words, the principal component is 
essentially constant with respect to sampling spin orien- 
tations. A small variance automatically implies approxi- 
mate conservation in the second sense but not necessarily 
in the first one. 

We emphasize that the interest here is not only in 
those PCs that have the smallest variances (and thus 
identify semi-conserved quantities) but also in those with 
the largest variances, which encode the most information 
about the inspiral dynamics. 



IV. VARIABLES FOR PCA 

In this paper we focus on analyzing the statistical prop- 
erties of differences between final and initial values of the 



following scalar products between the unit orbital angu- 
lar momentum and spin vectors in the following combi- 



nations: 



A(Si-L) 
A(S 2 -L) 

A(Si.Sa) 

A[(Sr • L)(S 2 • L)] 
A[(Sr.L) 2 ] 

A[(S 2 -L) 2 ] 



Si-L^-Si-Lli^AiL 

S 2 ■ L|/ — S 2 • L|i := A 2 l 

Si • S 2 |/ — Si ■ S 2 |i := A12 

(Si • t)\ f (S 2 ■ t)\ t - (Si • L)U(S 2 • t% := Ai L2L 

(§ 1 -t) 2 |/-(§i-l.) 2 |i:=A 1 i ilL 

(S 2 -L) 2 | / -(S 2 -L) 2 | i :=A 2L2L . 



(4.1) 
(4.2) 
(4.3) 
(4.4) 
(4.5) 
(4.6) 



r 



Other combinations are possible and might actually pro- 
vide further insight. 

In our numerical simulations and analytical calcula- 
tions we choose the initial spin orientations with uniform 
and uncorrelated probability distributions. The above 
scalar products are then also initially uncorrelated and 
their expectation values of (4.1)-(4.4) vanish, 



(Si-S 2 )L = <S 1 -L)| i = <S 2 -L)U = 



while those of (4.5 1 and (4.6 1 are 



((Si-IO^H&.L) 



1/3 



The orbital angular momentum and spin orientations 
naturally become correlated due to spin-orbit and spin- 
spin interactions as each of these binary black hole con- 
figurations evolve in time. However, at least within the 
PN approximation here considered, the orbital angular 
momentum and spin vectors remain perfectly uniformly 
distributed [5]. For example, a Kolmogorov-Smirnov test 
for a representative configuration returns a p-value of 
~ 10~ 5 when testing for lack of uniformness [8]. Higher 
PN expansions might introduce small biases [5] but if so 
they appear to be at a level in which approximating the 
mean of the above scalar products at any instant of time 
by zero is a very good approximation. 



V. INSTANTANEOUS APPROXIMATION 

Computing the covariance matrix requires sampling a 
large part of the space of initial spin orientations for each 
fixed set of binary black hole system parameters (tuples 
of masses and spin magnitudes) and necessitates comput- 
ing the solutions of the PN equations repeatedly for each 
set. Over long time periods (large tt — tj) we solve these 
equations numerically but we can also gain some insight 
into the structure of principal components by studying 
the evolution of the system over short durations. 

The changes in the scalar products used in our analysis 
can be calculated in the approximation where tf = ti+At 



for a small interval of time At. Using the equations of 
motion for the spin and orbital angular momentum and 
keeping terms through O(At) gives 



A 



1L 



3ra 2 x 2 w 2 A£ 



L-(Si xS 2 ) 



A 2L 
Ai 2 



^1L2L 



2M 

M- (Mw) 1/3 miXiSi 
Ail with lo2 

^Msixs.) 

2JV/4/3 v l > 

+ (Af W ) 1 / 3 (™ 2 xiSi-m 2 X2 S 2 

-^L.( Sl xS 2 ) 



Ai 



L1L 



miXiSi • L -m 2 x 2 S 2 • L 
3m 2 x 2 w 2 Ai ~ 



(Si • L) [L • (Si x S 2 ) 
m x xi(Ma;) 1/3 Si • L 



^2L2L 



M 
M 

Ailil with 1h2 



(5.1) 
(5.2) 

(5.3) 

(5.4) 

(5.5) 
(5.6) 



where the right sides are evaluated at the initial time 
(frequency) . 

In calculating the covariance of the scalar products we 
must first determine the expectation values of (5.1 1-(5.6) 



through 0{At). Doing this exactly is challenging be- 
cause At is found by solving the equations of motion for 
ui(t) and depend on the spin scalar products, which are 
stochastic variables, 



At = At(w,Si,S 2 ,L). 



(5.7) 



We perform this calculation because we use a common fi- 
nal frequency (as opposed to a common final time) in our 
analysis to provide a less gauge-dependent stopping crite- 
ria dependence. Physically, the interactions between the 
spins can cause a slight repulsion or attraction between 



the masses that increases or decreases, respectively, the 
elapsed time of the inspiral from w, to cuf. Therefore, 
we cannot simply factor At outside of the average nor 
are we able to solve the equations of motion analytically 
to find an expression for At. However, the spins affect 
the orbital frequency beginning at 1.5PN order and are 
small perturbations. Therefore, we can assume At is in- 
dependent of spin to a good approximation, and so the 
expectation values of the scalar products in (5.1 )-(5.6 ) 



all vanish in this approximation since the spins are un- 
corrected and have zero mean at the initial time. 

The elements of the covariance matrix are then cal- 



culated by forming the appropriate covariances of (5.1)- 
(5.6), which are all 0(At 2 ) to leading order. Since the 



spin directions are uniformly distributed over a 2-sphere 
then the expectation value of an odd number of spins 
vanishes and 



(si 1 



ST 



i 



AW 



) (5.8) 



(2n + l)!! 

for an even number of spins where the last ■ • • indicates 
all possible pairings of indices and k = 1,2. After some 
algebra we find that the elements of the covariance matrix 
are 



Ctj = (At) 2 C„ + <3(At 3 



(5.9) 



where the C^- are in some cases too lengthy to display in 
full detail. Since Cy = 0(At 2 ) at leading order then so 
also are the eigenvalues so that 

(5.10) 



CIV, = A , v . 



where A = (At) 2 A. We order the eigenvectors V\ and 
associated Ai such that Ai > A2 > ... > A n for an n x n 
covariance matrix. 



VI. SPIN-ORBIT VARIABLES (SO) 

We illustrate this approach in detail with a simple case 
and make contact with previous conservation results. We 
start building towards our more general result by first 
doing a PC A using only the two spin-orbit variables in 
p~Ij) and (|42l, 



A 1L := A(Si • L) , A 2L :=A(S 2 -L). 



(6.1) 



We remind the reader that we include spin-spin inter- 
actions in our simulations and analytical calculations by 
using the PN equations of motion described in Section [TT] 
and the Appendix, but in this section we only use the 



variables in (6.1 ) for the PCA 



For each black hole mass and spin magnitudes (mj , Xj ) 
the covariance matrix for the variables (|6.1[) is 



f Cov(A 1L ,A 1L ) Cov(A 1L ,A 2L ) 
- I Cov(A 2L ,A 1L ) Cov(A 2L ,A 2L ) 



(6.2) 



where the entries can come either from our numerical 
simulations or the instantaneous approximation. Wc 
then diagonalize C to find the principal components. 
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FIG. 1: A graphical representation of the principal compo- 
nents for the spin-orbit variables of Section |VI| and the nu- 
merical data of Ail and A2L for a binary black hole system 
with mi — 0.4 and maximal spin magnitudes. We show only 
10 points here but the PCs were computed with 10 random 
spin orientations. 



A. Numerical simulations 

From our numerical simulations we find that, sampling 
across many random initial spin orientations, each of the 
principal components has zero mean over time (to nu- 
merical accuracy), (A£?°) = 0, a consequence of the 
spin orientation distributions remaining highly uniform 
during the inspiral. 

Wc find A 2 to be in the range ~ 10 -9 to 10 -4 for the pa- 
rameters we sampled (toi e [0.1,0.9] and xi,2 € [0.1, 1]). 
Furthermore, A 2 grows with both spin magnitudes, which 
is expected, but increases as the equal mass case is ap- 
proached (we explain the reason for this behavior below) ; 
see Fig. [5] 

As an example with (wi,TO2,Xi,X2) = 

(0.4,0.6,1.0,1.0), Figure [I] shows a graphical repre- 
sentation of the principal components overlaid on a 
scatter plot of the Ail and A 2 l data from 1,000 of our 
100,000 numerical simulations using random initial spin 
orientations. Notice that the first PC, which points 
along the direction of the eigenvector Vi with the largest 
eigenvalue Ai, captures the largest variation in the data 
while the second PC, pointing along V 2 , indicates that 
there is very little spread in the data in that direction, 
which is also implied by the smallness of A 2 relative to 
Ai. 

The coefficients of the PCs, which are just the compo- 
nents of the corresponding eigenvectors, are functions of 
the black hole masses and spin magnitudes but also of 



their initial and final orbital frequencies, 

V,- = "Vj(mi,Xi>X2,Ui,Wf) , J = 1,2 . 

However, we have numerically found that the dependence 
on the initial and final frequencies is rather weak. This, 
together with the fact that A£ 2 has zero mean, implies 
that the function 

£f ° = t/ 1 §! • L + V£ S 2 ■ L 

is not only constrained to be nearly constant, in a sta- 
tistical sense, for all spin orientations for any black hole 
masses and spin magnitudes but also to be approximately 
constant as a function of time. This motivates and jus- 
tifies our instantaneous approximation described in Sec- 
tion M 



B. Instantaneous approximation 

In the instantaneous approximation we use the expres- 
sions in (5.1| and (5.2 1 to calculate the covariance ma- 



trix in (6.2 1 to leading order in At to find that Ci 



(At) 2 Cy + 0(At 3 ) with 



C12 = - ^m 1 m 2 XiX2^ i = C21 
C22 = C11 with 1 «-> 2. 



(6.3) 

(6.4) 
(6.5) 



We then diagonalize C to find the PCs as with our nu- 
merical simulations. We find that 



A£f° 
A£|° 



/ii(m 2 X2Ai L - miXiA 2L ) , (6.6) 

^ 2 A(S -L)=: A 0L , (6.7) 



with nx = [(mixi) 2 
Here, 



(«i2X2) 2 ] 1/2 and ^ 2 = m/M. 



More recently, it has also been shown that at the same 
PN order the quantity So • L/||L|| 2 is also a constant 
of motion if including quadrupole-monopole interactions, 



i.e. 



A(S -L/||L|| 2 )=0, (6.8) 

and ignoring spin-spin interactions and radiation reaction 



From (6.6) and (6.7) the instantaneous approximation 
implies that the PC with the smallest variance when in- 
cluding spin-spin and radiative effects in the equations 
of motion is neither A(So ■ L/||L|| 2 ) nor A e ffL but Aol 
instead [cf, Eqs. (6.7)]. The differences betwe en A l and 
A(S • L/||L|| 2 ) are detailed in Section | VIC 

Furthermore, the instantaneous eigenvalues (vari- 
ances) in the instantaneous approximation are given by 

A so_ 



(Ai) 2 A=° (j' = 1,2) with 



A* u = 



lOAf 4 / 3 



2 2, 2 2\ 

m x Xi +m 2 x 2 ) 



[5Af 4 / 3 ( 

+ ml m\ xt xl u 2/3 



;so 



™i ml xl X2 u 
IOM 4 / 3 



,14/3 



(6.9) 
(6.10) 



These expressions provide insight regarding the depen- 
dence of the variances on the black hole masses and 



spin magnitudes. For example, by analyzing Eq. (6.10) 



one can see that it increases with both spin magnitudes 
(which is expected) but also that it increases as the equal 
mass case is approached, as we already found from our 
numerical simulations (see Fig. [5| . 



C. Comparisons with other quantities 

AS f ° is the PC describing the largest variations in the 
data and A£f ° is essentially conserved, as both a func- 
tion of time and with respect to initial spin orientation 
variations, 



m 2 



nil 
m 2 



is one of the two spin vectors introduced in llj in the 
context of the effective one-body (EOB) formalism. In 
that reference it was also shown that at 2PN order in 
spin effects the scalar product between the effective spin 



Seff = I 1 



3m 2 
Ami 



Si 



1 



3m! 
Amy 



and the unit orbital angular momentum S e ff • L is a con- 
stant of motion if ignoring spin-spin interactions and ra- 
diative effects (see also [5]). That is, for each initial con- 
figuration this scalar product is constant in time, 

A cffL = A(S cff • L) = 0. 



cSO 
Co 



fi 2 Sq • L w constant. 



(6.11) 



In the approximations of Ref. [H] A£f° an d (6.8) are 
essentially equivalent to each other. Here, however, we 
are including radiative effects and spin-spin interactions 
and find that A£f i s conserved much better than So • 
L/||L|| 2 . 

For example, Figure [2] shows probability distributions 
from our numerical simulations for maximally spinning 
black holes with mi = 0.4. Displayed are the two numer- 
ical PCs A£f°,A£f° (solid lines), along with their in- 



stantaneous expressions (6.6), (6.7) (dotted lines), as well 
as the quantity in expression ( |6.8[ ) (dashed line). The nu- 
merical values and instantaneous approximations for the 
PCs are essentially indistinguishable. The variances for 
the PCs are Ai = 3.5 x 10~ 2 and A 2 = 3.9 x 10~ 5 , respec- 
tively, while the latter is to be compared with a variance 
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FIG. 2: Probability distributions for A8f° (wide) and AS- 
(narrow) for maximally spinning black holes with mi = 0.4. 
The solid lines are the numerical principal components and 
the dotted lines are principal components from the instanta- 



neous approximation in ( 6.6 \ and ( 6.7 1; note that the approxi 



mations are virtually indistinguishable from the numerical re- 
sul ts. T he dashed line shows the distribution of the quantity 
in (6.8 1. All the distributions have been normalized with re- 



spect to their maximum possible values so that the maximum 
range in these plots is [—1,1]. Notice that A£i° accounts for 
most of the variation in the data and AS f is sharply peaked 
around zero. 



of 6.9 x 10" 2 for (6.8) 



Notice that A£f ° has a consid- 



erably sharper distribution compared to expression (6.8) 
in Figure pi In fact, the variance of So • L/||L|| 2 is of 
the order of magnitude of the first principal component 
A£f°, which is the one describing the largest variation 



in the data. 

In Figure[2]we have not shown the probability distribu- 
tion for A c ffL because it is practically indistinguishable 
from A l- This is expected because So and S e ff are de- 
fined as sums of the individual spins with similar order 
of magnitude coefficients (in fact, So and S c g are propor- 
tional to each other in the equal mass and spin case). 



VII. 



SPIN-ORBIT AND SPIN-SPIN VARIABLES 

(SOSS) 



Section VI ignored spin-spin variables in the analysis. 
Here we extend the previous analysis using Ail and A2L 
by including A12 in the PCA. There are now three prin- 
cipal components, which are linear combinations of these 
variables, and as before we rank them by decreasing vari- 
ances. 

From our numerical simulations we find the smallest 
variance A3 has the same range as the smallest variance 



in Section VI However, the behavior is different with re- 
spect to the mass ratio. Namely, for any pair of spin 
magnitudes the variance A3 is smallest for small mass 
ratios, grows as the masses get comparable — taking a 
maximum around mi ~ 0.4 — and then decreases as the 
equal mass case is approached (see Fig. [5] for the maxi- 
mally spinning case xi — X2 — !)■ As in Section VI for a 



fixed mass ratio the variance is monotonically increasing 
with respect to both spin magnitudes. 

In the instantaneous approximation the infinitesimal 
covariance matrix is C = (At) 2 C + 0(At 3 ). with 



lOAf 6 



xS 2 



(5m 2 + xSl/M 2 ) 
5xS\S% 



T] 1 M 2 

X ^ S ^(5m 1 5m + xSl/M 2 ) 
r\ ' 



>1 



xS 2 
r] 2 M 2 
x l ' 2 S 1 



5xSiS 2 
V 
(5m 2 + xSl/M 2 ) 

(5m 2 Sm - xS 2 /M 2 



■''' S2 (bm^m + xSl/M 2 ) 



T) 

x 1 / 2 S 1 



(5m 2 6m - xS 2 /M 2 ) 



>1 

5M 2 Sm 2 +x(S 2 + S 2 ) 



(7.1) 



where Si = |S,| = m 2 \i, f] = mim^/Af 2 is the symmetric 
mass ratio, M — mi+m 2 is the total mass, Sm = m\—m 2 
and x = (Moj) 2/3 . 

The principal components and variances are straight- 
forward to compute but their expressions are rather 
lengthy so we only display those corresponding to the 
smallest eigenvalue. In this case 



\SOSS 
A 3 



with A| oss 



= (At) 2 \f oss 
0, and its associated (non- normalized) 



eigenvector is 



V 



SOSS 



(mlxi, m 2X2,m 1 xim 2 X2(Muj) 1/3 ) . (7.2) 



The third principal component is thus 
1 f(Muj) 1 / 3 



A£l oss 



with 



| V SOSS| 



77117712 



Si 



A(Si-S 2 ) + A(S- 



10 J , 

(7.3) 



the total spin angular momentum. In the instantaneous 
approximation the variance of A£f is zero (as is its 
expectation value) and equals the sum of non-negative 
numbers. It therefore follows that A£ | oss itself vanishes 



and that £ 



SOSSl _ £SOSS 



O(At), or 



(Mw) 1/3 A(Si • S 2 ) = -mim 2 A(S • L) + 0{At). (7.4) 

Similarly, in the extreme mass ratio limit, mi -C m 2 , 
(7.3) becomes 



A£ soss = !^l Xl ( m2W )i/3 A (s 1 . g 2 ) + A(S 2 • L) + O(At) 
m 2 

to leading order in the mass ratio. The fact that in our 
numerical simulations we find A£^ to be nearly conserved 
for finite times implies that 



(Afw) 1/3 A(Si • S 2 ) ~ -mim 2 A(S • L). 



VIII. SPIN-ORBIT AND SPIN-SPIN 
VARIABLES - NONLINEAR (NL) 



(7.5) 



Here we consider the six variables from Eqs. (4.1) 



(4.6), which include nonlinear dependence on the spin- 
orbit variables. When performing a PCA in these vari- 
ables we find that there are three principal components 
(A£f L , A£^ L and A£g h in order of decreasing vari- 
ance) with small variances. From our numerical simu- 
lations we find that these eigenvalues are all bounded 
from above by small numbers for any binary black hole 
mass and spin magnitudes. More precisely, for the span 
of frequencies considered in this paper they are in the 

1(T 9 to 1(T 4 and 



ranges A^ L - 1(T 7 to 1(T 2 

\NL 
A 6 



\NL 
A 5 



-7 

ict 11 to 10- 5 

As an example, Figure Wl shows the variance of A£ , 



NL 
6 

(i.e., the eigenvalue A^ L ) as a function of the black hole 
masses and spin magnitudes for evolutions up to a final 
frequency of w/ = 0.05. For any pair of spin magnitudes 
the variance is smallest for very different binary masses, 
grows as the latter get comparable, and then rapidly de- 
creases to zero as the equal mass case is approached. 
For a fixed mass ratio the variance is monotonically in- 
creasing with respect to both spin magnitudes. From 
Fig. [4] one can also notice that we can place the bound 
A^ L ^2x 10 -5 for any possible configuration of binary 
black hole masses and spin magnitudes. 



We emphasize that the PCs £^ L , £^ L and S, 



g L are dis- 



tinctly new quantities different from either So • L, S G ff • L 
or So • L/||L|| 2 . This can be noticed, for example, from 
Table [jl where we explicitly list the components of £^ L , 
£^ and £§ h for a particular mass and spin magnitude 
configuration. One can see that the spin-spin and the 
non-linear spin-orbit contributions are non-negligible. In 
addition, recall that the principal components are by con- 
struction statistically uncorrelated. 

The PCs £^ L , £^ L and £^ h are semi-conserved quan- 
tities in the two senses used in this paper in the presence 



TABLE I: Components of the normalized eigenvectors for the 
principal components with the three smallest variances for a 
binary black hole system with m\ — 0.4 (mi + rri2 = 1) and 
maximal spin magnitudes. 
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A 2L 


-0.408 


0.121 
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-0.604 


0.012 
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Ai_L2L 


-0.160 


0.757 


-0.132 


Aitiz, 


0.054 


0.248 


-0.018 


A2L2L 


0.224 


0.590 


-0.040 



of both radiative corrections (in the form of radiation 
reaction driving the inspiral) and spin-spin interactions. 
This is in contrast to S e ff • L and So ■ L/||L|| 2 , which 
are truly conserved when ignoring radiative effects and 
spin-spin interactions. However, in the presence of these 
corrections A ffL, A(So • L/||L|| 2 ) and Aol have larger 
variances than both A£f L and A£^ h . This is illustrated 
in Figure k| where we show the distributions of A£f L , 
A£|^ L and^f!^, normalized to their maximum possible 
value, for a binary black hole system with m\ — 0.4 and 
maximal spin magnitudes. This choice of parameters cor- 
responds to one of the largest variances for these three 
principal components and is thus indicative of a black 
hole binary with some of the "worst" statistics. We also 
show the normalized distributions for Aol and A e g l , 
which would be identically zero if ignoring spin-spin and 
radiative effects. However, they are not when includ- 
ing those effects and instead we find that while Aol and 
A c ff l are clustered around zero with a small variance of 
sa 5 x 10~ 5 the quantities A£^ L and A£^ L both have a 
variance smaller by a factor of « 10, indicating that £f h 
and £^ h are better conserved throughout the duration of 
the inspirals. 

Figure [5] shows the improvement obtained when in- 
cluding spin-spin variables in the PCA and when further 
including nonlinear spin-orbit ones. Shown is the small- 
est variance, for maximally spinning black holes, as a 
function of one of the black hole masses for the three PC 
analyses considered in this paper: 1) SO - purely spin- 
orbit terms, 2) SOSS - spin-orbit and spin-spin terms, 
and 3) Nonlinear - spin-orbit, spin-spin and nonlinear 
terms. The variance becomes smaller with each term in- 
cluded in the PCA. For the purely spin-orbit variables 
the smallest variance in the equal-mass case is near max- 
imum. In the other cases the equal mass configuration, 
in contrast, has the minimum variance with A^ L < 10~ 12 
for the PCA including the nonlinear variables. 



A. Binaries with equal masses and spin magnitudes 

Including the nonlinear spin-orbit variables into the 
principal component analysis reveals an interesting rela- 
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FIG. 3: Probability distribution of the three principal compo- 
nents with the smallest variance (A£f L (green), A£f h (blue) 
and A£^ L (red)). For comparison, Aol (solid black) and 
A e fl l (dashed), which are truly zero when ignoring spin-spin 
and radiative corrections, are also plotted. 
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tion in the specific case of a binary black hole with equal 
masses and spin magnitudes undergoing a quasi-circular 
inspiral. It turns out that, as hinted above, for this sub- 
set of parameter values A^ L is zero within numerical pre- 
cision, which indicates that not only does A£q L have 
a vanishing variance but, together with the observation 
that within such precision (A£ 6 ) is also zero, that A£§ L 
is itself zero. In other words, £g h is truly conserved, at 
least for the quasi-circular inspiral of binary black holes 
with equal masses and spin magnitudes. 

From the instantaneous approximation one can actu- 
ally see that £$ L is given by 



, NL _ 2S 1 -S 2 + (S 1 -L)(S 2 -L) 
' 6 " V5 



(8.1) 



and by using the PN equations of motion one can show 
that th e tim e derivative of £§ h vanishes. The conserva- 
tion of I 



1 1 implies the angle between Si and S 2 is fixed 



at every time by the angles of the spin orientations with 
the orbital angular momentum. 



B. Predictions using principal components 



Given that £|* L , 



pNL 
c 5 



and £g h 



are semi-conserved to 
varying degrees depending on the choice of binary black 
hole parameters it is tempting to use these quantities to 
predict the three relevant angles between the spin ori- 
entations and the orbital angular momentum. Figure [6] 
shows the errors in predicting Ail, A 2 l and A i2 for one 
of the cases with largest variances. The errors are de- 
fined as the predicted values minus those obtained from 
our numerical simulations. In calculating the predicted 
ones, we assume that A£f§ 6 are exactly zero and solve 



FIG. 4: Variance of the principal component AS g as a func- 
tion of the black hole spin magnitudes xi, X2 an d their masses 
(recall that throughout this paper the total mass is set to 
mi + 77i2 — 1 in all numerical simulations). 



for the quantity of interest in terms of the other vari- 
ables, the latter being replaced by their values from our 
numerical simulations. 

Figure shows that the assumption A£g L = does 
quite well at predicting either Ail or A 2 l and gives 
better predictions than assuming that either A l or 
A c ff l are zero. In predicting Ai 2 we see that assuming 
A£f L — gives the best predictions, albeit not partic- 
ularly good ones. Errors in approximating a principal 
component by exactly zero propagate in different ways 
into the scalar products when predicting the latter, ac- 
cording to their weights in the principal components (this 
is standard propagation of errors analysis). In particu- 
lar, the large errors in predicting Ai 2 are a consequence 
of the spin-spin interaction (and thus the corresponding 
weight in the principal component) being far weaker than 
spin-orbit corrections to the inspiral dynamics. 

The predictions in Figure [6] assume that one knows all 
variables but one with certainty. However, since we have 
three independent semi-conserved quantities, 



A£ 4 NL « 0, A£ 5 NL 



0, 



A£ r ^ L 



0. 



and three angles uniquely describing the orientations of 
the spin and orbital angular momenta vectors, it is in 
principle possible to predict these angles based solely on 
their distributions at the initial time. Unfortunately, in 
general this gives very large errors. For example, in the 
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FIG. 5: Smallest eigenvalue as a function of one of the 
black hole masses (maximally spinning, with total mass set 
to mi + vri2 = 1)- Shown are the purely spin-orbit (SO), the 
spin-orbit and spin-spin (SOSS), and the spin-orbit, spin-spin 
and nonlinear terms (nonlinear) cases. These correspond to 
eigenvalues A2 , A§ oss and Ag* 1 ", respectively. The equal- mass 
configuration has a near-maximum value for the SO case and 
a minimum in the other ones (with Ag < 10~ 12 ). 
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FIG. 6: Errors in the predictions of Ail, A2L and A12 for 
binary black holes with mi = 0.4 (total mass M — 1) and 
maximal spin magnitudes. The errors are defined as the pre- 
dicted values (found by solving A£f L = (green), A£f h = 
(blue) and A£q — (red)) minus those from our numerical 
simulations. For comparison, we have also plotted the errors 
in these predictions when imposing A(So-L) = (solid black) 
and A(Sbh ■ L) = (dashed). Note that none of these last 
two conditions can be used to make a prediction for A12. 



maximal spinning and equal mass case the errors in this 
prediction for Ail, A2L and A12 have, as probability 
distributions, standard deviations (not variances) of 0.32, 
0.32 and 0.05, respectively; and they are considerably 
larger for unequal masses. 



IX. GENERAL RELATIVISTIC DYNAMICS 

Since the PN approximation is highly accurate in the 
inspiral regime, one might expect that in the fully Gen- 
eral Rclativistic (GR) case the principal components with 
smallest variances discussed in this paper could also have 
small variances when computed with full GR data. If so, 



FIG. 7: Variances of A£ 2 so , A£f oss and A£g h as a function 
of time, using data from fully nonlinear GR simulations of 
binary black hole inspirals. 



a natural question would be whether that is still the case 
in the plunge regime. 

To have a preliminary sense of to what extent that 
might happen, we present three-dimensional numerical 
simulations of binary black holes using the full Einstein 
equations, starting around one and a half orbits before 
merger. Sampling the full parameter space of the binary 
black hole problem through numerical evolutions of the 
full Einstein equations is not possible. Even sampling a 
single configuration of masses and spin magnitudes would 
be unfeasible if 40,000 inspirals were needed. However, 
the fact that in the PN approximation the distributions 
of our principal components with smallest variances are 
so sharply peaked suggests that a small sample might be 
enough to reproduce their main features. For example, 
in the equal mass and equal spin case one can coarsen 
the sampling within the PN approximation to ~ 20 spin 
orientations and still reproduce the first significant digit 
of A£f oss when doing a PCA. 

We therefore present GR simulations of binary black 
hole configurations in an initial quasi-circular orbit with 
equal masses mi = mi = 0.5 and spin magnitudes 
Xi = Xi = 0.6, an initial separation of d — 6.2 M and 
20 random initial spin directions (Si, S2). These config- 
urations are evolved using the moving punctures tech- 
nique [T31 [H] all the way through merger. The spin vec- 
tors are measured over time following the isolated hori- 
zon approach [T5HT7] . The statistical results shown below 
are insensitive from removing some simulations from the 
analysis, hinting that such a small sample number might 
be sufficient to describe the salient features, just as in 
the PN approximation. 



In Figure \7\ we show the variances for AS. 



and AS, 



NL 
< 



SO 
2 ) 



A£f oss 



as defined by our PN PCA analysis, 



Eqs. (6.7), (|7.3[) and (8.1), but where the data for the 



scalar products (4.1 )-( 4.6 ) are taken from the fully rela- 



tivistic simulations. The variances for AS |° and AS§ oss 
appear indistinguishable from each other on the scale 
shown. The reason for this is that their difference is 



10 



proportional to the spin-spin interaction, which is small 
(see Eqs. (6.7), ( |7.3[ ) and Figure [5]). All the variances 
are considerably larger than the PN corresponding ones 
but the main features of the PN analysis seem to hold. 
The variances of Aff an d A£f oss are larger than the 
variance of AS$ h , as expected from our PC A analysis 
(see Fig. [5]), and its conservation properties (cf. the dis- 
cussion near Eq. (8.1)) at the PN level. In particular, 



notice that AS^ h is approximately constant (after the 
initial transients have settled for t <^ 25M) from about 
25M to 75M after which it grows exponentially during 
the plunge phase from about 75M to the final time shown 
in the figure. However, most importantly, the variance 
of our principal component AS^ h remains fairly small 
all the way until a common apparent horizon is found 
(which is just after the final time shown in Fig. uh. 

The results of these simulations suggest that to some 
extent the statistical structures identified in this paper 
within the PN approximation carry over to the inspiral 
and plunge regimes in the fully General Relativistic case. 



CONCLUSIONS 



three statistically uncorrelated semi-conserved quantities 
A£4,5,6 that do not correspond to any previously known 
expressions that we are aware of. These principal com- 
ponents are largely conserved when including spin-spin 
interactions and radiative corrections in the evolution to 
the extent that they are even reasonably conserved when 
using data from fully relativistic, three-dimensional sim- 
ulations of binary black holes in the inspiral and plunge 
regime, just up to merger. 

We expect these results to be useful in astrophysical 
simulations of binary hole inspirals and in generating 
waveform templates for gravitational wave data analy- 
sis. In particular, simplifying template generation may 
be realized with our PCA method by identifying those 
combinations of spin variables that are not only semi- 
conserved, and therefore effectively reduce the dimen- 
sionality of the parameter space (which may supplement 
the single-spin approximation used in [T8T 20 ) . but also 
those that encode the largest features in the data, and 
therefore are useful for efficiently sampling the parame- 
ter space thereby reducing the cost of constructing tem- 
plates. 



In summary, we have identified statistical constraints 
in binary black hole dynamics, from inspiral to merger, 
suggesting a dimensional reduction of the parameter 
space when probabilistic predictions are sufficient. We 
have also shown that these probabilistic predictions can 
be surprisingly accurate using numerical simulations or, 
to a lesser degree, using the analytical instantaneous ap- 
proximation. 

When using purely spin-orbit variables for a principal 
component analysis we have found that in the instan- 
taneous approximation the PC with smallest variance, 
A£.f (which is proportional to A(So • L)), is a statis- 
tically conserved quantity both in a statistical sense in 
terms of initial spin orientations but also as a function 
of time. We have also found the dominant component 
AS 1 ° describing the largest variation in the data. This 
demonstrates the ability of PCA to robustly identify con- 
served quantities as well as the most relevant features in 
the data. 

In our more general analysis we have identified 
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Appendix A: Post-Newtonian equations of motion 

The Post-Newtonian equations used in this paper are 
those of Ref. [§], Q33]), 



2 96 ,„, ^L 743 + 924?? , „ r n2/ o 



11 



^2 ( Xi t n -S i ( 11 ^ + 7r >l! )) -4, ! >/.: 



12 

i=l,2 

/ 34103 13661 



M 2 



59 



18 



V 18144 2016 



•-^-(4159 + 15876r])TT(Muj) 5/3 



r)+-ri z \ (M W ) 4 ' J - — »7XiX2 247(S! • S 2 ) - 721(L„ • Si)(L„ • S 2 ) (M« 



48 
/ 16447322263 1712 



n4/3 



16 



V 139708800 



105 



IE ■ 



^ +(- 



273811877 451 , 88- , 

1 it 0ri)n 

1088640 48 3 " ' 



541 2 5605 3 



-:n 



2592 



896 

j _ {MuY^dS 

n ~~ nM 2 ~dt 



856, ,„„,„,. ,o/o,\ /w x9 , 4415 358675 91495 2 , ,,, n7/ , 
; Zo3(16(Mw) 2/3 ) (Mw) 2 + ( h T] H 7? 2 )7r(Mu;) 7/3 



105 



4032 



6048 



1512 



(Al) 

(A2) 
(A3) 



where dS/dt = dSi/dt + dS 2 /dt, j E = 0.577. . . is Euler's can be computed via ||L„|| = ryM 5 / 3 ^ 1 / 3 . 

constant, and 9 — 1039/4620. The total mass is denoted 

by M = mi + fni and n — mim2/M 2 is the symmetric The evolution of the individual spin vectors Si for the 

mass ratio. The magnitude of the angular momentum 2 black holes is described by a precession around fij with 



J7i 
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and fl 2 is obtained by 1 <-> 2. 
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